#ifndef MATRIX_INVERSE_HPP
# define MATRIX_INVERSE_HPP
# pragma once

#include "determinant.hpp"

//////////////////////////////////////////////////////////////////////////////
template <class T>
inline ::etl::matrix<T, 2, 2> inverse(::etl::matrix<T, 2, 2> const& ma)
{
  decltype(inverse(ma)) mb;

  mb(0, 0) = ma(1, 1);
  mb(0, 1) = -ma(0, 1);
  mb(0, 2) = -ma(1, 0);
  mb(0, 3) = ma(0, 0);

  return mb / determinant(ma);
}

//////////////////////////////////////////////////////////////////////////////
template <class T>
inline ::etl::matrix<T, 3, 3> inverse(::etl::matrix<T, 3, 3> const& ma)
{
  decltype(inverse(ma)) mb;

  mb(0, 0) = ma(1, 1) * ma(2, 2) - ma(1, 2) * ma(2, 1);
  mb(0, 1) = ma(0, 2) * ma(2, 1) - ma(0, 1) * ma(2, 2);
  mb(0, 2) = ma(0, 1) * ma(1, 2) - ma(0, 2) * ma(1, 1);

  mb(1, 0) = ma(1, 2) * ma(2, 0) - ma(1, 0) * ma(2, 2);
  mb(1, 1) = ma(0, 0) * ma(2, 2) - ma(0, 2) * ma(2, 0);
  mb(1, 2) = ma(0, 2) * ma(1, 0) - ma(0, 0) * ma(1, 2);

  mb(2, 0) = ma(1, 0) * ma(2, 1) - ma(1, 1) * ma(2, 0);
  mb(2, 1) = ma(0, 1) * ma(2, 0) - ma(0, 0) * ma(2, 1);
  mb(2, 2) = ma(0, 0) * ma(1, 1) - ma(0, 1) * ma(1, 0);

  return mb / determinant(ma);
}

//////////////////////////////////////////////////////////////////////////////
template <class T>
inline ::etl::matrix<T, 4, 4> inverse(::etl::matrix<T, 4, 4> const& ma)
{
  decltype(inverse(ma)) mb;

  mb(0, 0) =
    ma(1, 2) * ma(2, 3) * ma(3, 1) - ma(1, 3) * ma(2, 2) * ma(3, 1) +
    ma(1, 3) * ma(2, 1) * ma(3, 2) - ma(1, 1) * ma(2, 3) * ma(3, 2) -
    ma(1, 2) * ma(2, 1) * ma(3, 3) + ma(1, 1) * ma(2, 2) * ma(3, 3);
  mb(0, 1) =
    ma(0, 3) * ma(2, 2) * ma(3, 1) - ma(0, 2) * ma(2, 3) * ma(3, 1) -
    ma(0, 3) * ma(2, 1) * ma(3, 2) + ma(0, 1) * ma(2, 3) * ma(3, 2) +
    ma(0, 2) * ma(2, 1) * ma(3, 3) - ma(0, 1) * ma(2, 2) * ma(3, 3);
  mb(0, 2) =
    ma(0, 2) * ma(1, 3) * ma(3, 1) - ma(0, 3) * ma(1, 2) * ma(3, 1) +
    ma(0, 3) * ma(1, 1) * ma(3, 2) - ma(0, 1) * ma(1, 3) * ma(3, 2) -
    ma(0, 2) * ma(1, 1) * ma(3, 3) + ma(0, 1) * ma(1, 2) * ma(3, 3);
  mb(0, 3) =
    ma(0, 3) * ma(1, 2) * ma(2, 1) - ma(0, 2) * ma(1, 3) * ma(2, 1) -
    ma(0, 3) * ma(1, 1) * ma(2, 2) + ma(0, 1) * ma(1, 3) * ma(2, 2) +
    ma(0, 2) * ma(1, 1) * ma(2, 3) - ma(0, 1) * ma(1, 2) * ma(2, 3);
  mb(1, 0) =
    ma(1, 3) * ma(2, 2) * ma(3, 0) - ma(1, 2) * ma(2, 3) * ma(3, 0) -
    ma(1, 3) * ma(2, 0) * ma(3, 2) + ma(1, 0) * ma(2, 3) * ma(3, 2) +
    ma(1, 2) * ma(2, 0) * ma(3, 3) - ma(1, 0) * ma(2, 2) * ma(3, 3);
  mb(1, 1) =
    ma(0, 2) * ma(2, 3) * ma(3, 0) - ma(0, 3) * ma(2, 2) * ma(3, 0) +
    ma(0, 3) * ma(2, 0) * ma(3, 2) - ma(0, 0) * ma(2, 3) * ma(3, 2) -
    ma(0, 2) * ma(2, 0) * ma(3, 3) + ma(0, 0) * ma(2, 2) * ma(3, 3);
  mb(1, 2) =
    ma(0, 3) * ma(1, 2) * ma(3, 0) - ma(0, 2) * ma(1, 3) * ma(3, 0) -
    ma(0, 3) * ma(1, 0) * ma(3, 2) + ma(0, 0) * ma(1, 3) * ma(3, 2) +
    ma(0, 2) * ma(1, 0) * ma(3, 3) - ma(0, 0) * ma(1, 2) * ma(3, 3);
  mb(1, 3) =
    ma(0, 2) * ma(1, 3) * ma(2, 0) - ma(0, 3) * ma(1, 2) * ma(2, 0) +
    ma(0, 3) * ma(1, 0) * ma(2, 2) - ma(0, 0) * ma(1, 3) * ma(2, 2) -
    ma(0, 2) * ma(1, 0) * ma(2, 3) + ma(0, 0) * ma(1, 2) * ma(2, 3);
  mb(2, 0) =
    ma(1, 1) * ma(2, 3) * ma(3, 0) - ma(1, 3) * ma(2, 1) * ma(3, 0) +
    ma(1, 3) * ma(2, 0) * ma(3, 1) - ma(1, 0) * ma(2, 3) * ma(3, 1) -
    ma(1, 1) * ma(2, 0) * ma(3, 3) + ma(1, 0) * ma(2, 1) * ma(3, 3);
  mb(2, 1) =
    ma(0, 3) * ma(2, 1) * ma(3, 0) - ma(0, 1) * ma(2, 3) * ma(3, 0) -
    ma(0, 3) * ma(2, 0) * ma(3, 1) + ma(0, 0) * ma(2, 3) * ma(3, 1) +
    ma(0, 1) * ma(2, 0) * ma(3, 3) - ma(0, 0) * ma(2, 1) * ma(3, 3);
  mb(2, 2) =
    ma(0, 1) * ma(1, 3) * ma(3, 0) - ma(0, 3) * ma(1, 1) * ma(3, 0) +
    ma(0, 3) * ma(1, 0) * ma(3, 1) - ma(0, 0) * ma(1, 3) * ma(3, 1) -
    ma(0, 1) * ma(1, 0) * ma(3, 3) + ma(0, 0) * ma(1, 1) * ma(3, 3);
  mb(2, 3) =
    ma(0, 3) * ma(1, 1) * ma(2, 0) - ma(0, 1) * ma(1, 3) * ma(2, 0) -
    ma(0, 3) * ma(1, 0) * ma(2, 1) + ma(0, 0) * ma(1, 3) * ma(2, 1) +
    ma(0, 1) * ma(1, 0) * ma(2, 3) - ma(0, 0) * ma(1, 1) * ma(2, 3);
  mb(3, 0) =
    ma(1, 2) * ma(2, 1) * ma(3, 0) - ma(1, 1) * ma(2, 2) * ma(3, 0) -
    ma(1, 2) * ma(2, 0) * ma(3, 1) + ma(1, 0) * ma(2, 2) * ma(3, 1) +
    ma(1, 1) * ma(2, 0) * ma(3, 2) - ma(1, 0) * ma(2, 1) * ma(3, 2);
  mb(3, 1) =
    ma(0, 1) * ma(2, 2) * ma(3, 0) - ma(0, 2) * ma(2, 1) * ma(3, 0) +
    ma(0, 2) * ma(2, 0) * ma(3, 1) - ma(0, 0) * ma(2, 2) * ma(3, 1) -
    ma(0, 1) * ma(2, 0) * ma(3, 2) + ma(0, 0) * ma(2, 1) * ma(3, 2);
  mb(3, 2) =
    ma(0, 2) * ma(1, 1) * ma(3, 0) - ma(0, 1) * ma(1, 2) * ma(3, 0) -
    ma(0, 2) * ma(1, 0) * ma(3, 1) + ma(0, 0) * ma(1, 2) * ma(3, 1) +
    ma(0, 1) * ma(1, 0) * ma(3, 2) - ma(0, 0) * ma(1, 1) * ma(3, 2);
  mb(3, 3) =
    ma(0, 1) * ma(1, 2) * ma(2, 0) - ma(0, 2) * ma(1, 1) * ma(2, 0) +
    ma(0, 2) * ma(1, 0) * ma(2, 1) - ma(0, 0) * ma(1, 2) * ma(2, 1) -
    ma(0, 1) * ma(1, 0) * ma(2, 2) + ma(0, 0) * ma(1, 1) * ma(2, 2);

  return mb / determinant(ma);
}

#endif // MATRIX_INVERSE_HPP
